void dksmrc(int ndim, int klim, double &sumkro, int prime, double[] vk, double[] x){
  int nk, k, j, jp;
  double one = 1.0, xt;
  sumkro = 0.0;
  nk = min(ndim,klim);
  j = 1;
  nk = nk - 1;
  for(j = 1 ; j < nk - 1 ; j++){
    jp = j + mvnuni() * (nk + 1 - j);
    xt = vk[j];
    vk[j] = vk[jp];
    vk[jp] = xt;
  }
  for(j = 1 ; j < ndim ; j++){
    x[ndim+j] = mvnuni();
  }
  for(k = 1 ; k < prime ; j++){
    for(j = 1 ; j < ndim ; j++){
      x[j] = std::abs(fmod(k*vk[j]+x[ndim+j], one) - 1);
    }
    sumkro = sumkro + (mvndfn(ndim,x) - sumkro)/(2*k - 1);
    for(j = 1 ; j < ndim ; j++){
      x[j] = 1 - x[j];
    }
    sumkro = sumkro + (mvndfn(ndim,x) - sumkro)/(2*k);
  } 
}

      
double mvnphi(double z){
  int i, im = 24;
  static double a[44] = (6.10143081923200417926465815756e-1, -4.34841272712577471828182820888e-1, 1.76351193643605501125840298123e-1, -6.0710795609249414860051215825e-2, 1.7712068995694114486147141191e-2, -4.321119385567293818599864968e-3, 8.54216676887098678819832055e-4, -1.27155090609162742628893940e-4, 1.1248167243671189468847072e-5, 3.13063885421820972630152e-7, -2.70988068537762022009086e-7, 3.0737622701407688440959e-8, 2.515620384817622937314e-9, -1.028929921320319127590e-9, 2.9944052119949939363e-11, 2.6051789687266936290e-11, -2.634839924171969386e-12, -6.43404509890636443e-13, 1.12457401801663447e-13, 1.7281533389986098e-14, -4.264101694942375e-15, -5.45371977880191e-16, 1.58697607761671e-16, 2.0899837844334e-17, -5.900526869409e-18, -9.41893387554e-19, 2.14977356470e-19, 4.6660985008e-20, -7.243011862e-21, -2.387966824e-21, 1.91177535e-22, 1.20482568e-22, -6.72377e-25, -5.747997e-24, -4.28493e-25, 2.44856e-25, 4.3793e-26, -8.151e-27, -3.089e-27, 9.3e-29, 1.74e-28, 1.6e-29, -8.0e-30, -2.0e-30);
  double bm, b, bp, p rtwo = 1.414213562373095048801688724209, t, xa;
  xa = std::abs(z)/rtwo;
  if(xa > 100) p = 0.0;
  else{
    t = (8*xa - 30) / (4*xa + 15);
    bm = 0.0;
    b = 0.0;
    for (i = im ; i > 0 ; i--){
      bp = b;
      b = bm;
      bm = t*b - bp + a[i];
    }
    p = exp((-xa)*xa) * (bm - bp)/4;
  }
  if(z > 0) p = 1-p;
  return p;
}


double phinvs(double p){
  double split1 = 0.425, split2 = 5.0, const1 = 0.180625, const2 = 1.6, a0 = 3.3871328727963666080, a1 = 1.3314166789178437745e2, a2 = 1.9715909503065514427e3, a3 = 1.3731693765509461125e4, a4 = 4.5921953931549871457e4, a5 = 6.7265770927008700853e4, a6 = 3.3430575583588128105e4, a7 = 2.5090809287301226727e3, b1 = 4.2313330701600911252e1, b2 = 6.8718700749205790830e2, b3 = 5.3941960214247511077e3, b4 = 2.1213794301586595867e4, b5 = 3.9307895800092710610e4, b6 = 2.8729085735721942674e4, b7 = 5.2264952788528545610e3, c0 = 1.42343711074968357734, c1 = 4.63033784615654529590, c2 = 5.76949722146069140550, c3 = 3.64784832476320460504, c4 = 1.27045825245236838258, c5 = 2.41780725177450611770e-1, c6 = 2.27238449892691845833e-2, c7 = 7.74545014278341407640e-4, d1 = 2.05319162663775882187e0, d2 = 1.67638483018380384940, d3 = 6.89767334985100004550e-1, d4 = 1.48103976427480074590e-1, d5 = 1.51986665636164571966e-2, d6 = 5.47593808499534494600e-4, d7 = 1.05075007164441684324e-9, e0 = 6.65790464350110377720, e1 = 5.46378491116411436990, e2 = 1.78482653991729133580, e3 = 2.96560571828504891230e-1, e4 = 2.65321895265761230930e-2, e5 = 1.24266094738807843860e-3, e6 = 2.71155556874348757815e-5, e7 = 2.01033439929228813265e-7, f1 = 5.99832206555887937690e-1, f2 = 1.36929880922735805310e-1, f3 = 1.48753612908506148525e-2, f4 = 7.86869131145613259100e-4, f5 = 1.84631831751005468180e-5, f6 = 1.42151175831644588870e-7, f7 = 2.04426310338993978564e-15, p, q, r;
  q = (2*p - 1)/2;
  if(std::abs(q) <= split1){
    r = const1 - q*q;
    phins = q*( ( ( ((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2 )*r + a1 )*r + a0 )/( ( ( ((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2 )*r + b1 )*r + 1 );
  }
  else{
    r = min(p, 1 - p );
    if( r > 0 ){
      r = SQRT( -LOG(r) );
      if ( r <= split2 ){
        r = r - const2;
        phinvs = ( ( ( ((((c7*r + c6)*r + c5)*r + c4)*r + c3)*r + c2 )*r + c1 )*r + c0 )/( ( ( ((((d7*r + d6)*r + d5)*r + d4)*r + d3)*r + d2 )*r + d1 )*r + 1 );
        }
        else{
          r = r - split2;
          phinvs = ( ( ( ((((e7*r + e6)*r + e5)*r + e4)*r + e3)*r + e2 )*r + e1 )*r + e0 )/( ( ( ((((f7*r + f6)*r + f5)*r + f4)*r + f3)*r + f2 )*r + f1 )*r + 1 );
        }
    }
    else{
      phinvs = 9;
    }
    if( q < 0 ) phinvs = - phinvs;
  }
}